 
 
!
!
!
!
!
!****************************************************************
!*                                                              *
!
!
!
!
!
!dk geom
      subroutine gauss3v(npar, radius, prob, r, p, deriv2, tmp) 
!-----------------------------------------------
!   M o d u l e s 
!-----------------------------------------------
      USE vast_kind_param, ONLY:  double 
      use modify_com_M 
 
!...Translated by Pacific-Sierra Research 77to90  4.3E  14:13:36   8/20/02  
!...Switches: -yf -x1             
      implicit none
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      integer  :: npar 
      real(double)  :: radius(*) 
      real(double)  :: prob(*) 
      real(double)  :: r(*) 
      real(double)  :: p(*) 
      real(double)  :: deriv2(*) 
      real(double)  :: tmp(*) 
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      integer :: n, icount 
      real(double) :: pi, rmax, dr, p0, ptotl, derivl, derivr, dprob 
!-----------------------------------------------
!
!=======================================================================
!
!
!=======================================================================
!
      pi = acos(-1.) 
!
!-----------------------------------------------------------------------
!
!     compute probability of having particles inside a sphere of radius
!     r for a 1d3v Maxwellian distribution function.  Note that the
!     distribution function is cutoff beyond three standard deviations.
!
!      r and p are dummy arrays
!
!
      rmax = 3.0 
      dr = rmax/float(npar) 
      p0 = 4./sqrt(pi) 
      ptotl = 0.0 
      do n = 1, npar 
         r(n) = float(n)*dr 
         ptotl = ptotl + p0*r(n)**2*exp((-r(n)**2))*dr 
         p(n) = ptotl 
      end do 
      derivl = 1./(p0*r(1)**2*exp((-r(1)**2))) 
      derivr = 1./(p0*r(npar)**2*exp((-r(npar)**2))) 
      dprob = p(npar)/float(npar) 
      do n = 1, npar 
         prob(n) = float(n)*dprob 
      end do 
!
!-----------------------------------------------------------------------
!
!     From a tabulated array of probability as a function of speed, wher
!     speed is equally spaced between 0 and three standard deviations, w
!     use cubic spline to interpolate and to compute speed as a function
!     of probability, where probability is now equally spaced between 0
!     and 1.
!
      call spline (p, r, npar, derivl, derivr, deriv2, tmp) 
      call splint (p, r, deriv2, npar, prob, radius) 
      do n = 2, npar 
         radius(n) = 0.5*(radius(n)+radius(n-1)) 
      end do 
      radius(1) = radius(1)/2. 
      icount = 0 
      return  
      end subroutine gauss3v 
